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optimal  design  of  morphing  structures  comprised  of  active  smart  materials 
through  locally  controllable  actuation  and  activation.  The  computational  optimal 
design  framework  for  smart  material  morphing  structures  was  then  extended  to 
substantially  improve  the  computational  efficiency  of  the  overall  solution 
procedure  by  replacing  the  non-gradient-based  optimization  with  gradient-based 
optimization.  In  addition,  this  work  sought  to  explore  the  use  of  reduced-order 
modeling  for  accurate  and  efficient  inverse  problem  solution  algorithms,  and 
investigated  inverse  problem  solution  efficiency  in  a  general  sense.  The  reduced- 
order  modeling  component  was  primary  focused  on  nondestructive  evaluation- 
type  inverse  problems  (e.g.,  inverse  characterization).  A  computational  approach 
was  established  to  create  physics-based  reduced-order  models  for  representing  the 
behavior  of  a  boundary  value  problem  (e.g.,  solid  mechanics,  heat  transfer,  etc.)  in 
a  generally  applicable  way,  with  minimal  computational  cost,  but  with  accuracy 
commensurate  with  traditional  finite  element  analysis  methods.  Then,  an  entirely 
new  variation  on  the  use  of  proper-orthogonal  decomposition  (POD)  bases  for 
inverse  material  characterization  problems  through  a  Gappy  POD  reconstruction 
procedure  combined  with  direct  inversion  was  hypothesized,  developed,  and 
tested.  Lastly,  the  use  of  multi-objective  optimization,  rather  than  the  typically 
utilized  single-objective  format,  for  non- gradient-based  optimization-based 
inverse  problem  solution  strategies  was  investigated  and  evaluated.  Through 
simulated  examples,  the  multi-objective  approach  was  shown  to  maintain  a 
substantially  higher  diversity  in  the  potential  solution  set  throughout  the  non¬ 
gradient-based  optimization  in  comparison  to  single-objective  approaches. 
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Accomplishments/New  Findings: 


A  computational  framework  was  established  for  the  optimal  design  of  morphing 
structures  comprised  of  active  smart  materials  through  locally  controllable 
actuation  and  activation.  The  framework  combined  finite  element  analysis  to 
predict  the  morphing  process  for  a  given  set  of  activation  and  actuation 
parameters  with  nonlinear  optimization  to  identify  the  optimal  set  of  those 
parameters.  For  context  the  framework  was  developed  for  the  specific  application 
of  morphing  structures  comprised  of  thermally  responsive  polymers,  and  therefore 
the  activation  was  applied  thermal  stimulus.  The  framework  was  capable  of 
considering  linear  transient  heat  transfer  coupled  with  quasi-static  solid 
mechanics.  Specific  targets  of  the  design  process  that  could  be  considered 
including  the  desired  motion  or  kinematic  of  the  structure  as  well  as  the  energy 
and  time  required  to  complete  the  morphing  process. 

The  initial  incarnation  of  the  computational  optimal  design  framework  for  smart 
material  morphing  structures  utilized  non-gradient-based  optimization  to  ensure 
that  the  design  solution  was  near  to  the  global  optimum  regardless  of 
computational  expense.  This  framework  was  numerically  evaluated  through 
simulated  case  studies.  The  framework  was  shown  to  be  capable  of  identifying 
design  parameters,  including  activation  and  actuation  locations  and  sequencing,  to 
produce  substantial  savings  in  terms  of  total  energy  cost  and  more  accurately 
achieve  a  desired  shape  change  in  comparison  to  standard  intuitive  designs. 

The  computational  optimal  design  framework  for  smart  material  morphing 
structures  was  extended  to  substantially  improve  the  computational  efficiency  of 
the  overall  solution  procedure  by  replacing  the  non-gradient-based  optimization 
with  gradient-based  optimization.  Furthermore,  to  even  further  improve  the 
computational  efficiency  an  adjoint  approach  was  formulated  and  implemented  to 
calculate  the  gradient  of  the  design  objectives  with  respect  to  the  morphing 
system  activation  and  actuation  parameters  that  only  required  the  equivalent  of 
two  finite  element  analyses  per  iteration.  Again,  the  design  procedure  with 
gradient-based  optimization  was  shown  to  be  capable  of  identifying  design 
parameters,  including  activation  and  actuation  locations  and  sequencing,  to 
produce  substantial  savings  in  terms  of  total  energy  cost  and  more  accurately 
achieve  a  desired  shape  change  in  comparison  to  standard  intuitive  designs. 
Furthermore,  although  the  design  solutions  were  less  optimal  than  those  that  could 
be  obtained  with  non-gradient-based  optimization,  the  computational  cost  of  the 
gradient-based  approach  was  multiple  orders  of  magnitude  less  than  that  of  the 
gradient-based  optimization. 

A  computational  approach  was  established  to  create  physics-based  reduced-order 
models  (ROM)  for  representing  the  behavior  of  a  boundary  value  problem  (e.g., 
solid  mechanics,  heat  transfer,  etc.)  in  a  generally  applicable  way,  with  minimal 
computational  cost,  but  with  accuracy  commensurate  with  traditional  finite 
element  analysis  methods.  The  approach  utilized  the  reduced-basis  ROM 
technique  and  was  established  in  the  context  of  inverse  problem  applications,  such 
that  the  accuracy  of  the  ROM  is  maximized  over  a  range  of  system  parameters 
that  will  ultimately  be  solved  for  through  an  inverse  problem  solution  strategy. 

The  ROM  creation  approach  was  tested  with  example  nondestructive  evaluation 


problems  and  was  shown  to  provide  a  sufficiently  accurate  representation  of  the 
system  physics  to  be  utilized  in  an  inverse  problem  solution  technique,  and 
required  only  a  fraction  of  the  cost  of  a  standard  finite  element  analysis  method. 

An  entirely  new  variation  on  the  use  of  proper-orthogonal  decomposition  (POD) 
bases  for  inverse  material  characterization  problems  was  hypothesized,  developed, 
and  tested.  This  new  variation  creates  a  POD  basis  from  a  previously  generated 
or  observed  set  of  system  response  fields  for  the  system  to  be  characterized  as  is 
typical.  However,  rather  than  use  the  basis  to  recreate  a  reduced-order  model,  the 
basis  is  used  with  the  Gappy  POD  machine  learning  technique  to  create  a  tool  that 
predicts  a  full- field  response  from  partial  field  measurements  (as  would  be 
available  from  nondestructive  testing  techniques).  Then,  the  full-field  response  is 
utilized  with  a  direct  inversion  strategy  to  predict  the  inverse  characterization 
solution.  The  approach  was  shown  through  simulated  nondestructive  evaluation 
problems  to  require  a  fraction  of  the  cost  of  typical  iterative  inverse  problem 
solution  techniques,  to  maintain  general  applicable,  and  to  still  provide  accurate 
inverse  solution  estimates. 

The  use  of  multi-objective  optimization,  rather  than  the  typically  utilized  single¬ 
objective  format,  for  non-gradient-based  optimization-based  inverse  problem 
solution  strategies  was  investigated  and  evaluated.  Through  simulated  examples, 
the  multi-objective  approach  was  shown  to  maintain  a  substantially  higher 
diversity  in  the  potential  solution  set  throughout  the  non-gradient-based 
optimization  in  comparison  to  single-objective  approaches.  This  increase  in 
diversity  was  shown  to  substantially  improve  the  efficiency  of  the  optimization- 
based  solution  strategy  as  well  as  the  consistency  and  accuracy  of  the  final 
solution  estimate.  In  addition,  a  strategy  was  developed  to  utilize  this  solution 
diversity  throughout  the  inverse  problem  solution  procedure  to  adaptive  modify 
the  parameterization  of  the  unknown  field  of  the  inverse  problem,  thereby  further 
increasing  the  accuracy  and  efficiency  of  the  inverse  problem  solution  estimation 
procedure. 


SUMMARY 


A  strategy  to  obtain  maximally  efficient  and  accurate  morphing  structures  composed  of  ac¬ 
tive  materials  such  as  shape  memory  polymers  (SMP)  through  synchronization  of  adaptable  and 
localized  activation  and  actuation  was  established  and  numerically  tested.  A  computational  in¬ 
verse  mechanics  framework  was  created  that  combines  a  computational  representation  of  the  SMP 
thermo-mechanical  behavior  with  a  nonlinear  optimization  algorithm  to  determine  location,  mag¬ 
nitude  and  sequencing  of  the  activation  and  actuation  to  obtain  a  desired  shape  change  subject  to 
design  objectives  such  and  prevention  of  damage.  The  concept  of  localized  activation  along  with 
the  optimal  design  strategy  were  shown  to  be  able  to  produce  far  more  energy  efficient  morph¬ 
ing  structures  and  more  accurately  reach  the  desired  shape  change  in  comparison  to  traditional 
methods  that  require  complete  structural  activation  prior  to  actuation. 

The  computational  strategy  for  estimation  of  the  optimal  parameters  relating  to  the  distribution 
and  sequencing  of  activation  and  actuation  for  a  morphing  smart  material  structure  or  structural 
component  to  efficiently  and  effectively  achieve  a  desired  morphing  function  was  extended  to  sub¬ 
stantially  increase  computational  efficiency.  In  particular,  a  formulation  of  the  adjoint  method  was 
developed  to  be  utilized  to  estimate  the  gradient  of  the  objective  function(s)  with  respect  to  the 
activation  and  actuation  parameters  with  minimal  computational  expense  within  a  gradient-based 
optimization  strategy,  which  then  provides  an  optimization  process  that  is  substantially  compu¬ 
tationally  efficient  overall.  Overall,  the  computational  optimal  design  approach  with  the  adjoint 
method  was  shown  to  provide  the  capability  to  efficiently  identify  activation  and  actuation  pa¬ 
rameters  to  achieve  desired  morphing  capabilities.  Furthermore,  the  computational  approach  was 
shown  to  be  capable  of  determining  energy-efficient  design  solutions  for  a  diverse  set  of  target  shape 
changes  with  fixed  instrumentation,  providing  the  potential  for  substantial  functionality  beyond 
what  could  be  expected  through  traditional  empirical  design  strategies. 

A  generally  applicable  algorithm  for  adaptive  generation  of  data  ensembles  to  efficiently  cre¬ 
ate  accurate  computational  mechanics  reduced-order  models  (ROM)  for  use  in  computational  ap¬ 
proaches  to  approximate  inverse  problem  solutions  was  developed  and  numerically  evaluated.  The 
ROM  approach  considered  was  based  on  identifying  the  optimal  low- dimensional  basis  to  be  used 
within  a  Galerkin  weak-form  finite  element  method  to  provide  substantially  reduced  computational 
cost  while  maintaining  accuracy  relative  to  that  of  a  (traditional)  full-order  finite  element  model. 
The  core  hypothesis  of  the  algorithm  presented  is  that  maximizing  the  diversity,  as  defined  in  a 
measurable  sense,  of  the  full-order  models  used  to  create  the  ROM  will  maximize  the  accuracy  of 
the  ROM  over  a  range  of  input  system  parameters.  The  adaptive  snapshot  generation  algorithm 
was  shown  to  produce  ROMs  that  can  accurately  estimate  response  fields  over  a  wide  range  of 
input  parameters,  and  which  are  substantially  more  accurate  than  ROMs  created  from  randomly 
generated  snapshot  sets.  Moreover,  the  accurate  generalization  of  the  adaptively  generated  ROMs 
was  shown  to  be  sufficient  to  consistently  produce  accurate  inverse  characterization  solution  esti¬ 
mates  with  a  fraction  of  the  computational  expense  that  would  be  required  to  do  so  with  full-order 
analyses. 

A  new  and  unique  approach  for  computationally  efficient  inverse  material  characterization  from 
partial-field  response  measurements  that  combines  the  Gappy  proper  orthogonal  decomposition 
(POD)  machine  learning  technique  with  a  physics-based  direct  inversion  strategy  was  created 
and  evaluated.  Gappy  POD  is  used  to  derive  a  data  reconstruction  tool  from  a  set  of  potential 


system  response  fields  that  are  generated  from  available  a  priori  information  regarding  the  potential 
distribution  of  the  unknown  material  properties.  Then,  the  Gappy  POD  technique  is  applied 
to  reconstruct  the  full  spatial  distribution  of  the  system  response  from  whatever  portion  of  the 
response  field  has  been  measured  with  the  chosen  system  testing  method.  Lastly,  a  direct  inversion 
strategy  is  presented  that  is  derived  from  the  equations  governing  the  system  response  (i.e.,  physics 
of  the  system),  which  utilizes  the  full- field  response  reconstructed  by  Gappy  POD  to  produce  an 
estimate  of  the  spatial  distribution  of  the  unknown  material  properties.  The  inversion  procedure 
was  shown  to  have  the  capability  to  efficiently  provide  accurate  estimates  to  material  property 
distributions  from  partial-field  response  measurements.  The  direct  inversion  with  Gappy  POD 
response  estimation  was  also  shown  to  be  substantially  tolerant  to  noise  in  comparison  to  the 
direct  inversion  given  measured  full-field  response. 

A  multi-objective  optimization-based  computational  approach  to  nondestructive  evaluation  of 
material  properties  was  numerically  evaluated.  The  multi-objective  approach  provides  a  substan¬ 
tial  improvement  in  the  capabilities  to  traverse  the  optimization  search  space  to  minimize  the 
measurement  error  and  produce  accurate  damage  estimates.  In  addition,  a  novel  self-evolving  pa¬ 
rameterization  approach  for  nondestructive  evaluation  was  developed  and  numerically  evaluated. 
The  adaptive  approach  utilizes  the  substantial  solution  diversity  that  is  uniquely  provided  by  multi¬ 
objective  optimization  to  iteratively  build  up  the  parameterization  and  accurately  characterize  all 
property  changes  with  the  minimum  dimensional  parameterization.  The  nondestructive  evalua¬ 
tion  approach  with  self-evolving  parameterization  was  shown  to  provide  an  accurate  and  efficient 
process  for  the  solution  of  inverse  characterization  problems. 


Optimal  Design  of  Locally  Activated  Smart  Material  Morphing  Structures: 

Figure  1  illustrates  the  general  form  of  the  morphing  structure  design  (and  analogously  control) 
problems  considered.  In  this  concept,  the  morphing  process  begins  with  some  portion  of  the  domain 
heated  with  a  controllable  transient  temperature  distribution  and/or  surface  heat  flux  ( TA(x,t ) 
or  qA(x,t)).  Once  a  sufficient  portion  of  the  structure  is  softened  through  activation,  mechanical 
actuation  begins  through  controlled  transient  displacement  and/or  force  ( uA(x,t )  or  rA(x,t ))  to 
deform  the  structure  into  a  desired  shape.  Therefore,  the  design/control  problem  to  identify  the 
activation  and  actuation  distributions  to  achieve  the  desired  shape  change  can  be  cast  in  a  general 
way  in  the  form  of  the  following  constrained  optimization  problem: 
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where  x  is  the  spatial  position  vector,  t  is  the  time,  Q  is  the  spatial  domain  of  the  structure,  tF 
is  the  total  time  of  the  morphing  process,  ||  •  ||nt  is  the  metric  norm  with  respect  to  the  spatial 
domain  and  time,  star9et  is  the  desired  shape  change,  sapprox  is  the  approximate  expected  shape 
change  given  the  design  parameters  predicted  by  a  numerical  representation  of  the  structure,  .JF 
and  J\y  are  the  total  thermal  energy  and  mechanical  work,  respectively,  predicted  by  the  numerical 
representation  to  perform  the  morph,  w  represents  additional  design  objectives  that  may  be  desired, 
cr  and  T  are  the  approximate  expected  internal  stress  tensor  and  internal  temperature,  respectively, 
predicted  by  the  numerical  representation  throughout  the  structure  during  the  morphing  process, 
cr(-)  is  a  suitable  scalar  stress  transform  (e.g.,  von  Mises  criterion),  and  adama9e  and  Tdama9e  are  the 
stress  and  temperature  limits,  respectively,  to  avoid  damaging  the  material.  A  framework  was  built 
to  solve  the  constrained  optimization  problem  (Eqns.  (1)  -  (2))  through  a  computational  inverse 
problem  solution  approach  that  combines  a  numerical  representation  of  the  morphing  system  (i.e. , 
forward  problem)  and  a  nonlinear  optimization  strategy  to  identify  the  parameters  that  minimize 
the  desired  objectives  while  satisfying  the  required  constraints. 


Figure  1:  Schematic  of  the  morphing  structure  design  problem  in  which  the  applied  activation  ( TA 
or  qA )  and  actuation  (ilA  or  ta)  are  to  be  determined  to  achieve  a  desired  shape  change. 


Two  optimization  approaches  were  investigated  for  the  solution  of  this  (now  PDE-constrained) 
optimization  problem.  The  first  portion  of  the  investigation  (i.e.,  the  proof-of-concept)  used  a 
genetic  algorithm,  a  stochastic  global  search  algorithm  that  mimics  the  process  of  evolution  in 
nature  [1,  2],  Genetic  algorithms  (GAs)  are  heuristic  optimization  approaches  that  generally  rely 
on  three  main  operations  to  evolve  a  population  of  potential  solutions  to  the  optimization  problem 
parameters:  survival,  reproduction,  and  mutation.  GAs  have  seen  substantial  use  in  recent  years 
due  to  their  ease  of  implementation  (e.g.,  there  is  no  need  for  complicated  gradient  calculations), 
global  search  capabilities,  and  simplicity  of  parallelization.  However,  GAs  may  not  be  ideal  for 
many  applications  in  computational  inverse  mechanics  due  to  the  large  number  of  function  evalu¬ 
ations  (i.e.,  numerical  simulations)  that  are  typically  required  to  converge  to  an  optimal  solution. 
To  create  a  substantially  more  realistically  applicable  framework,  the  second  portion  of  the  work 
utilized  gradient-based  optimization,  particularly  relying  on  the  adjoint  approach  to  expedite  gra¬ 
dient  calculations.  A  quasi-Newton  gradient-based  algorithm,  the  BFGS  interior  point  algorithm 
[3,  4,  5],  was  chosen  for  this  purpose. 

To  present  the  adjoint  formulation  developed  and  implemented,  consider  the  following  specific 
example  form  of  the  optimal  design/control  objective  functional: 

m  =  j(p)+c-  s(p>  (3) 

where  p  is  the  vector  of  unknown  design  parameters  to  be  determined  through  the  optimization 
(i.e.,  activation  and  actuation  parameters)  and  C  is  the  scalar  weighting  constant  (note  that  the 
specific  value  of  the  weighting  constant  allows  emphasis  of  the  optimization  effort  to  be  placed 
on  a  specific  objective,  energy  or  shape,  depending  on  the  relative  importance  for  a  particular 
application).  The  total  energy  of  the  morphing  process  (J  =  Jp  +  Jw)  can  be  defined  with: 


Jt  =  I  I  pc[T(x,  t )  —  To]<5(t  —  tp)dxdt 

It  Jo 


(4) 


and 

Jw  —  <r(if,  t)  ■  n(x )  •  u(x,  t)dxdt , 

where  8  is  the  Dirac  delta  function,  and  the  morphing  error  could  be  defined  as: 

M 


S  —  f  /  \u(x,t)  —  utar9et(x,  f)]2  S(x  —  x*m)8(t  —  t p)dxdt, 

•It  Jo.  m=1 


(5) 


(6) 


where  u  and  'utarflei  are  the  predicted  and  desired  (i.e.,  target)  displacement  fields,  respectively, 
{^mlm= i  is  *4ie  set  °f  M  discrete  spatial  locations  where  the  target  displacement  is  desired  to  be 
achieved  by  the  morph,  and  tp  is  the  time  at  the  completion  of  the  morphing  process.  For  simplicity 
of  illustration,  the  activation  process  can  be  assumed  to  occurred  entirely  through  a  temperature- 
dependent  Young’s,  and  the  behavior  of  the  structures  was  assumed  to  be  defined  by  linear  transient 
heat  transfer  and  linear  elastic  quasi-static  solid  mechanics.  Therefore,  the  morphing  process  of  a 


smart  material  structure  could  be  represented  by  the  following  initial  boundary  value  problem: 


kS72T(x,  t)  =  pc 
V  •  cr(f,  t )  =  0, 

T(x,  0)  =  T0, 

T(f,  t )  =  TA(x,  t ), 

—k\7T(x,  t )  •  n{x)  =  (p4, 

«(f,  t)  =  uA(x,  t), 
cr(f ,  t)  •  n(f)  =  rA(f,  A), 
cr(f,  A)  =  CIV  (E(T),v)  :  e(f,  t), 
e(f,  t)  =  |[Vu(x,  t)  +  V«(f,  t)T], 


Vfefi,  te  [o,tF], 

Vf  ell,  i  e  [o, 

V  x  G  id,  t  =  0, 

V  x  G  rT,  t  G  [0,  tF], 

V  x  G  r„  t  G  [0,tF], 

V  if  G  Tj/,  t  G  [0, £p], 
V  f  6  b,  16  [0,  tF], 
Vf  6  O,  t  G  [0,tF], 
Vf  6  fi,  16  [0,tF], 


(7) 


where  A:  is  the  thermal  conductivity,  p  is  the  mass  density,  c  is  the  specific  heat,  a  is  the  Cauchy 
stress  tensor,  e  is  the  small  strain  tensor,  CIV  is  the  fourth  order  elastic  stiffness  tensor,  E  is  the 
Young’s  modulus,  u  is  the  Poisson’s  ratio,  u  is  the  displacement  vector,  Q  is  the  spatial  domain,  T 
is  the  domain  boundary,  n  is  the  unit  outward  normal  vector  on  the  boundary  Fr  is  the  portion  of 
the  domain  boundary  where  the  temperature  is  specified,  Tq  is  the  portion  of  the  domain  boundary 
where  the  heat  flux  is  specified,  Tp  is  the  portion  of  the  domain  boundary  where  the  displacement 
is  specified,  and  is  the  portion  of  the  domain  boundary  where  the  traction  is  specified.  Following 
the  procedure  for  the  adjoint  method  [6],  for  the  above  forward  problem  the  adjoint  problem  for 
the  Lagrange  multipliers  (A  and  (p)  can  be  written  as: 


k  v2  X(x,t)  +  pc^§^-  -  !§f(v<F)  =  Vf  Gfi,  t  G  [0  ,tF], 

V  •  <rp(x,t)  +  CJ2m= 1 -  v?ar9et]S(x  -  x *m)S(t  -  tF)  =  0,  V  x  G  Q,  t  G  [0 ,tF], 


A(x,  tF)  —  0, 

V  f  G  12,  t  =  tp, 

A(f,  t)  —  0, 

V  i  6  Tj1,  t  G  [0,tF], 

—k  v  A(f,  t )  •  n(x)  =  0, 

V  x  6  Tq,  t  G  [0 ,tF], 

T 
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II 

v  x  g  rF,  t  G  [0,tF]. 

cTipix,  t )  •  n{x)  -  rA5(t  -  tF)  4  =  0, 

V  x  G  r?,  t  G  [0,  bp] 

Then,  the  gradient  of  the  objective  function  can  be  calculated  as: 

Ws  =  ftfrT (~k  V  A  •  n)  %^dxdt  -  ft  fVq  (X^)dxdt 
+  It  Jrv  (V  +  ™A)  i^dxdt  +  ft  fFu(-crp  ■  n  •  +  cra  •  n  •  )dxdt , 
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where 
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(V  '  «)I, 
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2(1  +  v)L  J  (l  +  u)(l-2 u) 

and  I  is  the  identity  tensor.  The  algorithm  for  the  calculation  of  the  gradient  using  the  adjoint 
method  can  be  summarized  as  follows: 


(1)  Solve  the  forward  initial  boundary  value  problem  (Eqns.  (7))  to  obtain  the  displacement  field 
u  and  the  temperature  field  T. 

(2)  Solve  the  adjoint  initial  boundary  value  problem  (Eqns.  (8))  to  obtain  the  Lagrange  multipliers 
A  and  <p. 

(3)  Substitute  u,  T,  A  and  (p  into  Eqn.  (9)  to  compute  the  gradient  of  the  objective  function. 

For  this  work,  the  forward  initial  boundary  value  problem  and  the  adjoint  problem  were  solved 
using  the  traditional  Galerkin  finite  element  method. 

To  test  the  inherent  benefits  and  challenges  of  the  developed  computational  approach,  numerical 
investigations  were  performed  with  example  design  applications  from  concepts  of  morphing  skeletal 
structural  components  (i.e.,  framing  and  connecting  elements).  A  smart  link  concept  with  direct 
shape  control  was  tested  to  verify  the  computational  framework,  as  well  as  to  compare  the  efficiency 
and  effectiveness  of  the  computational  approach  using  the  adjoint  method  with  the  non-gradient- 
based  genetic  algorithm  approach.  Then,  an  example  of  a  morphing  structural  backbone  with 
substantially  increased  complexity  of  design  parameters  and  objectives  was  considered.  As  shown 
in  Figure  2,  the  smart  link  consisted  of  a  20 mm  x  5 mm  x  1  mm  homogeneous  rectangular  prism 
composed  of  thermally  responsive  SMP.  The  morphing  procedure  considered  involved  first  heating 
the  structure  in  its  cast  configuration  to  initiate  the  activation  process.  After  some  time  of  purely 
heating,  actuation  would  be  applied  to  deform  the  structure  into  the  desired  shape.  As  such, 
the  total  time  to  perform  the  structural  morphing,  tp,  was  simply  the  summation  of  the  purely 
heating  time  prior  to  actuation,  tp,  and  the  actuation  time,  Im-  A  schematic  of  the  structural 
backbone  considered  for  the  second  example  is  shown  in  Figure  3.  The  structure  was  taken  to  be 
a  homogeneous  half-circular  prism  with  a  150mm  outer  radius  composed  of  thermally  responsive 
SMP.  The  morphing  procedure  applied  followed  a  similar  sequencing  pattern  as  the  first  example. 
However,  to  provide  greater  degrees  of  freedom  for  the  design/control  of  the  morphing  system, 
several  more  heating  pads  and  actuators  were  employed  in  this  case  compared  to  the  first  example, 
with  a  total  of  nine  pairs  of  approximately  45mm  long  heating  pads  and  eight  approximately 
15mm  long  actuation  regions,  all  equally  spaced  along  the  surface  of  the  structure.  Additionally, 
the  operation  of  the  structure  and  the  design/control  objectives  were  assumed  to  be  symmetric 
with  respect  to  the  horizontal  access.  Thus,  only  the  top  half  to  the  structure  needed  to  be  analyzed 
and  the  morphing  process  of  the  structure  was  defined  by  10  independent  heating  pads  and  four 


Figure  2:  Schematic  of  the  smart  link  design  concept. 


I 


Figure  3:  Schematic  of  the  Morphing  Structural  Backbone  Design  Concept. 

Table  1:  Control,  optimized  design  parameters  from  Genetic  Algorithm  optimization  (GA  Opt), 
optimized  design  parameters  from  quasi-Newton  gradient-based  optimization  using  the  adjoint 
method  (Adjoint  Opt),  resulting  energy  consumption  to  perform  the  morph,  and  the  percent  dif¬ 
ference  (%Diff)  with  respect  to  the  control  design  for  the  Smart  Link  example  ( mJ  =  10 ~3  Joule). 


Design 

Case 

t\ 

(s) 

(s) 

t'^P 

(s) 

tj* 

00 

Jt 

(mJ) 

Jw 

(mJ) 

J 

(mJ) 

Control 

80.0 

80.0 

80.0 

80.0 

14037.5 

10.4 

14047.9 

GA  Opt 

44.4 

42.5 

41.5 

69.3 

9231.2 

83.2 

9314.5 

%Diff 

-44.5 

-46.9 

-48.2 

-13.3 

-34.2 

+704.0 

-33.7 

Adjoint  Opt 

51.0 

57.3 

65.6 

66.7 

10833.9 

86.7 

10920.6 

%Diff 

-36.2 

-28.4 

-18.0 

-16.6 

-22.8 

+737.1 

-22.3 

actuators  (note  that  the  heating  pad  pairs  were  not  linked  in  this  example  as  was  the  case  in  the 
first  example). 

Table  1  shows  a  set  of  “control”  design  parameters,  the  optimized  design  parameters  from 
the  GA  optimization,  and  the  optimized  design  parameters  from  the  gradient-based  optimization 
using  the  adjoint  method,  along  with  the  corresponding  activation  thermal  energy  (Jr),  actuation 
mechanical  work  (Jw),  and  total  energy  (./),  and  the  percent  difference  for  these  energies  with 
respect  to  the  control  design  for  the  smart  link  example.  The  GA  results  represent  an  approximation 
of  the  best  possible  results  that  could  obtained  from  the  computational  design  procedure,  yielding 
an  energy  savings  of  more  than  30%  with  respect  to  the  control  design,  but  required  an  extensive 
computational  cost,  on  the  order  of  104  function  evaluations,  to  achieve  convergence.  Alternatively, 
the  gradient-based  optimization  results  were  able  to  yield  less  energy  savings  than  the  GA,  but  still 
a  substantial  amount  at  over  20%  with  respect  to  the  control  design,  but  required  two  orders  of 
magnitude  less  function  evaluations,  on  the  order  of  102,  than  the  GA  to  achieve  convergence  (even 


with  the  10  runs  with  different  initial  parameter  sets).  Furthermore,  it  is  often  the  case  for  design 
or  control  problems  such  as  this  that  an  approximation  of  the  global  minimum  is  unnecessary,  and  a 
local  solution  would  be  sufficient.  In  such  instances,  the  substantial  improvement  in  computational 
savings  provided  by  the  gradient-based  approach  with  the  adjoint  method  would  likely  outweigh  the 
loss  in  global  search  capabilities.  Table  2  shows  three  morphing  target  shapes  considered  (elliptic, 
square,  and  step-type),  the  optimal  design  solutions  for  each  target  shape,  and  the  corresponding 
temperature  distributions  at  the  initiation  of  morphing,  final  stress  distributions,  and  final  deformed 
shapes  for  the  morphing  structural  backbone  example.  All  of  the  optimal  design  solutions  can  be 
seen  to  have  produced  recognizable  approximations  to  the  target  shapes,  which  is  particularly 
significant  considering  that  each  of  the  three  shapes  is  considerably  different  from  one  another 
and  yet  the  morphing  mechanism  (i.e.,  activators  and  actuators)  used  to  achieve  the  shapes  is 
identical.  However,  there  was  a  degradation  in  the  morphing  accuracy  as  the  complexity  of  the 
target  shape  change  increased,  which  is  not  necessarily  unexpected.  This  morphing  accuracy  can  be 
seen  quantitatively  through  the  relative  difference  between  the  optimized  morphing  error  (Sopt)  and 
the  morphing  error  for  the  starting  shape  (So),  which  was  approximately  90%  for  the  elliptic  shape, 
40%  for  the  square  shape,  and  10%  for  the  step-type  shape.  The  target  step-type  shape  was  clearly 
the  most  challenging  to  achieve,  both  intuitively  and  based  on  the  optimization  results,  particularly 
since  it  was  the  only  shape  that  required  a  change  from  convex  to  concave.  Furthermore,  the  results 
would  imply  that  to  better  achieve  a  target  shape  similar  to  the  step-type  the  starting  shape  or  the 
morphing  mechanisms  (i.e.,  activators  and  actuators)  would  need  further  modification.  However, 
overall,  the  multiple  localized  activations  and  actuations  showed  the  capability  to  accurately  and 
efficiently  achieve  a  diverse  set  of  shape  changes  with  a  fixed  set  of  morphing  mechanisms,  and  the 
computational  approach  was  particularly  well-suited  to  facilitate  the  implementation,  particularly 
where  intuitive  design  approaches  would  be  infeasible. 


Table  2:  Optimized  (Opt)  design  solutions  (temperature  distribution  and  final  stress  distribution) 
with  respect  to  three  different  target  shapes  (ellipse,  square  and  step-type)  for  the  Morphing 
Structural  Backbone  example. 


Target 

Shape 

Temperature 
Distribution  (K) 

Stress 

Distribution  (Pa) 

S0  =0.380,  Sopt  =0.015 

l 

r  403 

t400 

\  .  375 

It 

'l  359 

6.4  x  104 

i4.0xl04 

=2.0  xlO4 

1  1.2  x  102 

S0  =0.109,  Sopt  =0.063 

n 

*  403 

#  K  400 

~  375 

1  E 

“  359 

_ 6.4  xlO4 

\  -;4.0xl04 

1  =  2.0xl04 

1  i 

1.2  X  102 

i 

|  S0 =1.055,  Sopt  =0.954 

ft 

...  403 

1  400 

4  ^375 

it 

■  359 

_  6.4  xlO4 

i4.0xl04 

\  =2.0  x  104 

1  i 

1.2  X  102 

Maximized  Reduced- Order  Model  Accuracy  for  Inverse  Problem  Applications: 

Although  generally  applicable,  the  following  discussion  of  an  approach  to  create  optimally 
accurate  physics-based  reduced-order  models  (ROM)  is  presented  within  the  context  of  steady-state 
harmonic  solid  mechanics  of  heterogeneous  solids  with  a  range  of  potential  material  parameters. 
Assuming  that  the  solid  considered  is  excited  harmonically  to  a  steady-state  in  the  linear  elastic 
range,  and  therefore  the  system  variables  vary  harmonically  in  time  with  angular  frequency  u>, 
and  neglecting  body  forces,  the  governing  equations  and  boundary  conditions  (i.e. ,  boundary  value 
problem)  from  conservation  of  momentum  can  be  written  as: 

V  •  cr(x,  u>)  +  u>2p(x)u(x,  u) )  =  0,  \/x  G  fl, 

cr(x,u> )  •  n(x )  =  T(x,u>),  \/x  G  IV, 

u(x,u>)  =  v?(x,u>),  \/x  G  Tf/,  (12) 

cr(x,cu )  =  CIV  :  e(T,  w), 

e(x,cu)  =  \  (V«(f,w)  +  Vu(x,co)T)  , 

where  x  is  the  spatial  position  vector,  cr(x,u)  is  the  stress  tensor,  p(x)  represents  the  density  of 
the  solid,  u(x,u> )  is  the  steady-state  displacement  amplitude  field,  T(x,u>)  is  the  applied  traction 
amplitude  vector,  m°(x,  lo)  is  the  applied  displacement  amplitude,  CIV  is  the  fourth-order  elasticity 
tensor,  Q  is  the  domain  of  the  solid,  n(x)  is  the  unit  outward  normal  vector  to  the  surface  of  the 
domain,  F,  and  Tt  and  r (}  are  the  portions  of  the  domain  surface  where  traction  and  displacement 
are  applied,  respectively,  such  that  Txljr^  =  F  and  Fr  f~)  F(/  =  0.  The  standard  weak  form 
Galerkin  approach  was  employed  for  this  work  to  approximate  the  solution  of  the  boundary  value 
problem  described  by  Eqns.  (12)  using  an  arbitrary  approximation  function  of  the  steady-state 
harmonic  displacement  field.  As  such,  the  weak  form  of  the  steady-state  dynamic  solid  mechanics 
problem  can  be  expressed  as: 

/  V5u(x)  :  cr(x,uj)dx  —  /  ui2p{x)5u{x)  •  u(x,  uj)dx  —  /  5u(x)  •  T(x,  uj)dx  =  0,  (13) 

Jq  Jq  JTt 

where  8u(x)  is  an  arbitrary  weight  function  vector  (i.e.,  trial  function  vector  or  virtual  displacement 
vector).  Therefore,  all  that  is  necessary  to  complete  the  Galerkin  approach  is  to  substitute  an 
approximation  for  the  displacement  field  and  the  weight  function  (using  the  same  basis  for  both) 
to  obtain  a  discretized  form  and  a  linear  system  of  equations  for  each  excitation  frequency.  The 
common  finite  element  approach  would  be  to  discretize  the  spatial  domain  into  elements  and  use 
polynomial  approximations  within  each  element  to  discretize  then  assemble  a  system  of  equations. 
However,  as  is  commonly  known,  this  finite  element  approach  (referred  to  as  the  full-order  modeling 
approach  herein)  typically  requires  at  least  many  thousands  of  degrees  of  freedom,  even  for  relatively 
simple  two-dimensional  problems  to  accurately  represent  the  physics  of  the  system.  Alternatively, 
the  objective  of  the  reduced-basis  form  of  reduced-order  modeling  is  to  identify  a  basis  that  is 
optimal  in  some  sense  for  representing  the  physics  of  the  system  under  consideration  with  far  fewer 
degrees  of  freedom  than  the  full-order  (i.e.,  traditional  finite  element)  model  (FOM). 

The  core  hypothesis  of  the  reduced-basis  reduced  order  modeling  approach  considered  in  the 
present  work  is  that  a  relatively  small  number  of  full-order  (i.e.,  traditional  finite  element)  analyses 
based  upon  different  values  of  the  input  parameters  of  interest  (material  parameters  herein)  contain 


fundamental  information  about  the  potential  solution  fields  of  the  BVP  and  can  be  used  to  derive 
a  low-dimensional  basis  that  can  predict  the  solution  fields  for  a  range  of  input  parameters  (not 
just  the  specific  parameter  values  used  to  generate  the  set  of  full-order  analyses)  with  reasonably 
sufficient  accuracy.  The  POD  approach  was  applied  to  derive  this  basis  from  a  set  of  full-order 
analyses.  POD  specifically  derives  the  low-dimensional  basis  such  that  the  difference  between  the 
original  full-order  data  and  the  best  approximation  to  that  data  with  this  basis  is  minimized  in 
an  L-2  average  sense  [7].  The  critical  question  is  then  how  to  select  the  set  of  input  parameters 
used  to  create  the  set  of  full-order  analyses  to  use  for  creating  the  POD  basis.  In  particular,  this 
dataset  must  be  generated  in  such  a  way  to  limit  the  number  of  full-order  simulations  necessary  to 
ensure  sufficiently  accurate  generalization  of  the  reduced-order  model  over  the  admissible  range  of 
the  input  parameters  of  interest.  The  problem  of  finding  a  set  of  snapshots  to  create  an  accurate 
ROM  can  be  outlined  in  a  general  way  as  follows.  Suppose  that  F  :  X  — »  V  is  a  finite  element  (i.e., 
full-order)  operator  that  maps  the  space  of  input  parameters  (e.g.,  material  properties)  X  to  the 
space  of  displacement  response  fields  V.  The  goal  is  to  find  a  subset  Xn  C  X  with  F  :  Xn  — >  Vn 
to  create  a  POD  reduced-order  model  operator  Vn)  :  X  — »  V,  such  that: 


if  F  :  x  —¥  V\  and  F^iVn)  :  x  — >•  v2,  then  ||ui  —  u2||  <  e,  \/x  £  X ,  and  Vvi,v2  £  V, 

where  e  is  some  acceptable  error  tolerance  for  the  ROM.  Extending  the  work  in  [7],  the  core 
hypothesis  of  the  proposed  adaptive  snapshot  generation  method  is  that  maximizing  the  diversity 
of  the  snapshots  created  within  the  space  of  the  input  parameters  of  interest  will  maximize  the 
generalization  of  the  resulting  reduced-order  model  over  that  parameter  space. 

Given  a  set  of  n0  input  parameters  sets  {7?}-=i  and  the  corresponding  full-order  analysis  re¬ 
sponse  fields  (i.e.,  snapshot),  a  measure  of  the  total  diversity  of  the  ith  snapshot  (created  with 
input  parameter  set  7*)  within  the  set  can  be  calculated  as: 

no 

=  (14) 

3= 1 


where 


j) 


«(®,7<)IU2(n)  •  ll«(£,7i)IUa(n)’ 


(15) 


Then,  provided  with  the  diversity  metric  for  each  snapshot  in  the  set  {R*  a  surrogate  model 

approach  can  create  an  approximate  mapping  between  the  input  parameters  and  the  diversity 
metric.  Any  preferred  machine  learning  technique  can  be  used  to  generate  the  surrogate  model, 
with  this  work  using  support  vector  regression  [8,  9].  The  surrogate  model  of  the  diversity  RSM( 7) 
can  then  be  easily  minimized  to  estimate  the  optimal  next  set  of  input  parameters  (within  the 
domain  of  the  parameters  X)  to  use  with  the  full-order  model  to  generate  another  set  of  snapshots 
that  would  maximize  the  diversity  of  the  snapshot  set  as: 


Minimize  RSM(  7). 


(16) 


Note  that  since  the  computational  cost  of  the  surrogate  model  is  negligible,  then  essentially  any 
preferred  global  optimization  algorithm  can  be  used,  regardless  of  the  algorithm  efficiency,  with  a 


Figure  4:  Flowchart  describing  the  adaptive  snapshot  generation  algorithm. 

genetic  algorithm  being  used  for  this  work.  Due  to  the  expectation  of  some  loss  of  accuracy  in  the 
surrogate  model,  an  additional  constraint  on  the  parameter  sets  was  added  to  the  surrogate  model 
optimization  for  the  present  work.  To  ensure  that  the  input  parameter  sets  are  not  overly  clustered 
in  the  parameter  space,  the  new  parameter  set  was  constrained  to  be  a  specified  minimum  distance 
5  from  every  other  parameter  set,  such  that: 

||7“7ill><*  for  i  =  1,2,...,  no,  (17) 

where  ||  •  ||  is  the  standard  /2-norm.  Once  the  new  set  of  input  parameters  is  determined  from  the 
solution  of  Eqns.  (16)  and  (17),  the  input  parameter  set  is  used  with  the  full-order  model  to  generate 
a  new  snapshot.  Figure  4  shows  a  flowchart  that  describes  the  overall  procedure  for  the  adaptive 
snapshot  generation  algorithm  to  maximize  snapshot  diversity  and  resulting  reduced-order  model 
accuracy.  At  the  completion  of  the  snapshot  generation  algorithm  a  POD  reduced-order  model  can 
be  generated  and  utilized. 

To  investigate  the  capabilities  and  potential  applicability  of  the  method  for  adaptive  genera¬ 
tion  of  optimal  reduced-order  models,  simulated  case  studies  were  considered  regarding  efficient 
and  accurate  modeling  of  the  deformation  of  structural  members  with  semi-localized  Young’s  mod¬ 
ulus  distributions.  The  examples  also  examined  the  capabilities  to  then  inversely  characterize  such 
material  property  distributions  using  a  computational  inverse  solution  procedure  relying  on  this 
modeling.  POD  ROMs  were  generated  through  the  adaptive  approach  to  maximize  diversity  of 
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Figure  5:  Schematic  for  NDE  plate  example  (note  that  the  ‘x’s  represent  the  sensor  locations). 


the  snapshot  sets,  and  the  accuracy  of  these  ROMs  was  quantified  with  respect  to  full-order  anal¬ 
ysis  (i.e.,  traditional  finite  element  analysis)  and  compared  to  the  accuracy  of  similar  POD  ROMs 
created  from  randomly  generated  snapshot  sets.  Then,  the  adaptively  generated  ROM  was  incor¬ 
porated  into  an  inverse  characterization  solution  procedure.  The  case  studies  considered  structures 
tested  using  frequency-response-based  nondestructive  testing  (NDT)  to  then  determine  the  mate¬ 
rial  properties  with  nondestructive  evaluation  (NDE).  The  NDT  consisted  of  a  localized  harmonic 
actuation  applied  normal  to  the  surface  of  the  structure  at  a  given  excitation  frequency  and  the 
resulting  steady-state  harmonic  vertical  displacement  amplitude  was  measured  at  a  set  of  discrete 
sensor  locations.  One  example,  shown  schematically  in  Figure  5,  consisted  of  a  lm  x  lm  x  0.02m 
aluminum  plate  subject  to  a  lk  Pa  harmonic  load  applied  to  a  5 cm  region  normal  to  the  top  surface 
of  the  plate,  excited  to  steady-state  with  an  actuation  frequency  of  400Hz.  The  material  behavior 
was  assumed  to  be  linear  elastic  with  a  homogeneous  density  and  Poison’s  ratio  of  2700  kg/m 3  and 
0.3,  respectively,  and  the  semi-localized  Young’s  modulus  distribution  was  assumed  to  be  defined 
with  a  radial  basis  function  (RBF)  as: 


E(x)  =  E0 


1  —  a  •  exp 


z-CII 


2 


5 


(18) 


where  E0  is  the  base  Young’s  modulus,  a  is  the  percentage  of  the  reduction  in  Young’s  modulus,  (' 
is  the  center  of  the  RBF,  and  c  is  the  breadth  of  the  RBF.  In  other  words,  each  Young’s  modulus 
distribution  considered  was  parametrized  by  four  parameters,  such  that  7  =  [«,  £1,  £2,  c]T.  The 
base  Young’s  modulus  was  assumed  to  be  known  as  a  standard  nominal  value  for  aluminum,  such 
that  E0  =  69 GPa. 

For  the  adaptive  snapshot  generation  approach,  an  initial  random  set  of  10  snapshots  was 
generated,  and  then  the  adaptive  surrogate  modeling  approach  was  iteratively  applied  to  generate 
the  remaining  snapshots  in  the  set  used  to  create  the  ROM.  To  examine  the  dependence  on  the 
total  number  of  snapshots,  snapshot  sets  of  10  (i.e.,  the  original  randomly  generated  set),  20,  30, 


40,  50,  and  60  were  investigated,  in  turn.  The  forward  modeling  accuracy  of  the  ROMs  was  first 
tested  directly  in  comparison  to  the  full-order  modeling  for  100  randomly  generated  parameter  sets 
that  were  not  included  in  the  snapshot  sets,  to  directly  quantify  the  generalization  capabilities  of 
the  ROMs  before  considering  an  inverse  characterization  problem.  To  test  the  ROM  accuracy  a 
standard  relative  error  metric  was  utilized  for  each  parameter  set  as  follows: 


Error (7) 


\u 


\ ROM 


\\uFOM(x,uj,^)\\a 


(19) 


where  uROM  and  uFOM  are  the  displacement  response  fields  calculated  with  the  ROM  and  full-order 
model,  respectively,  and  ||  •  ||q  is  the  chosen  norm  over  the  spatial  domain,  O,  with  both  the  L2 
and  Loo  norms  being  considered  in  the  following.  To  provide  a  baseline  for  comparing  the  accuracy 
of  the  adaptively  generated  ROMs,  ROMs  were  also  created  for  the  examples  using  an  equivalent 
total  number  of  randomly  generated  snapshots  (generated  in  the  same  format  as  the  initial  set  for 
the  adaptive  approach).  To  test  the  efficacy  of  the  resulting  ROMs  to  be  used  in  a  computational 
inverse  problem  solution  procedure,  each  example  case  considered  a  corresponding  set  of  tests 
in  which  the  ROMs  were  used  in  an  NDE  procedure  to  estimate  the  stiffness  parameters  of  the 
structures  given  simulated  NDT  measurements  for  several  test  cases.  The  inverse  problem  was 
cast  as  an  optimization  problem  to  determine  the  material  parameters  that  minimize  the  relative 
difference  between  the  simulated  experimental  NDT  measurements  and  the  response  predicted  by 
the  ROM  as: 


Minimize 

'yGX 


X7l  {uROM(Xi,UJ,  7)  -  Uexp(Xi,Uj)Y 
£?=1 


(20) 


where  again  X  is  the  domain  of  the  unknown  stiffness  parameters  and  uexp  is  the  simulated  experi¬ 
mental  displacement  responses.  A  standard  genetic  algorithm  was  again  applied  to  solve  the  above 
optimization  problem  and  identify  the  parameters  to  estimate  the  Young’s  modulus  distributions, 
and  therefore,  estimate  the  solution  to  the  inverse  problem.  The  quality  of  the  final  inverse  prob¬ 
lem  solution  estimates  were  quantified  through  the  relative  L2-error  between  the  Young’s  modulus 
distribution  defined  by  the  parameters  used  to  create  the  simulated  experimental  data  and  that 
estimated  by  the  inverse  characterization  results  as: 


^  fQ  (E(x,  7exp)  —  E(x,  7 inv))2  dx^j 
(fn(E(xXexp'l)2dx)1/2 


(21) 


where  yexp  are  the  parameters  used  to  create  the  simulated  experimental  measurement  data  and 
7 mv  are  the  corresponding  inverse  solution  estimates. 

Figure  6  shows  the  average  and  standard  deviation  of  the  relative  ROM  error  for  the  100  test 
cases  for  both  the  adaptively  generated  ROMs  and  the  randomly  generated  ROMs.  As  would  be 
expected,  for  both  approaches,  the  average  error  as  well  as  the  standard  deviation  of  the  error  for 
the  resulting  ROM  decreased  as  the  number  of  snapshots  used  to  construct  the  ROM  increased. 
More  interestingly,  the  ROM  error  corresponding  to  the  adaptively  generated  snapshots  was  sub¬ 
stantially  lower  than  the  the  ROM  error  corresponding  to  the  randomly  generated  snapshots  by 
approximately  a  factor  of  2  or  more  for  every  size  of  the  snapshot  set.  In  addition,  the  standard 
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Figure  6:  Average  and  standard  deviation  (error  bars)  with  respect  to  the  100  test  cases  of  the 
relative  L2  and  ROM  errors  for  the  randomly  generated  (Random)  and  the  adaptively  generated 
(Adaptive)  ROMs  for  the  plate  example. 


deviation  of  the  error  for  the  adaptively  generated  ROMs  decreased  considerably  more  quickly 
than  the  randomly  generated  ROMs,  and  while  the  adaptively  generated  ROMs  appeared  to  have 
converged  to  some  degree  in  terms  of  the  error  at  50  snapshots,  the  randomly  generated  ROMs 
show  no  such  signs  of  convergence.  Table  3  shows  the  material  parameters  used  to  create  the 
simulated  experimental  measurement  data  for  the  five  scenarios  considered  and  the  corresponding 
parameters  estimated  by  the  inverse  solution  process  with  the  various  ROMs  (built  from  20,  40, 
and  60  total  snapshots),  the  respective  ROM  measurement  error  for  the  inverse  solutions  (as  de¬ 
fined  by  Eqn.  (20)),  the  respective  FOM  measurement  error  for  the  inverse  solutions  (substituting 
the  FOM  in  place  of  the  ROM  in  Eqn.  (20)),  and  the  error  in  the  respective  Young’s  modulus 
distributions  predicted  by  the  inverse  solution  estimates  (as  defined  by  Eqn.  (21)).  In  addition,  to 
provide  further  perspective  on  the  relative  accuracy  of  the  Young’s  modulus  distributions  obtained 
by  the  inverse  solution  process,  Figure  7  shows  (as  a  representative  example)  the  target  (i.e.,  sim¬ 
ulated  experimental)  Young’s  modulus  distribution  for  a  material  property  scenario  compared  to 
the  Young’s  modulus  distribution  that  was  inversely  estimated  using  the  ROM  built  from  60  adap¬ 
tively  generated  snapshots.  Overall,  the  optimization  process  was  able  to  sufficiently  match  the 
ROM  response  to  the  measurement  data.  More  importantly,  the  FOM  responses  with  the  inverse 
solution  estimates  also  sufficiently  matched  the  measurement  data,  even  though  the  optimization 
was  performed  with  the  ROM.  Moreover,  the  FOM  measurement  error  was  minimally  higher  than 
the  ROM  measurement  error.  In  other  words,  the  inverse  problem  solution  estimates  obtained 
with  the  ROMs  were  nearly  as  accurate  with  respect  to  the  FOM  in  terms  of  the  measurement 
data,  and  were  still  within  an  error  range  in  terms  of  the  FOM  to  be  considered  an  inverse  problem 
solution  estimate.  As  would  be  expected,  corresponding  to  the  accuracy  in  the  measurement  error, 
the  resulting  estimates  of  the  Young’s  modulus  distributions  were  accurate  for  all  three  ROMs, 
with  Young’s  modulus  reconstruction  errors  of  less  than  1%  for  every  test. 


Table  3:  Target  (i.e.,  simulated  experimental)  values  for  the  RBF  amplitude  (a),  the  breadth  of 
the  RBF  (ci),  and  the  horizontal  and  vertical  locations  of  the  center  of  the  RBF  (Ci,C2)  defining 
the  Young’s  modulus  distribution,  the  corresponding  parameters  estimated  with  the  inverse  char¬ 
acterization  process  using  the  ROMs  created  with  20  (ROM-20),  40  (ROM-40),  and  60  (ROM-60) 
adaptively  generated  snapshots,  and  the  respective  ROM  measurement  error  (Me),  the  FOM  mea¬ 
surement  error  (Fe),  and  the  error  in  the  predicted  Young’s  modulus  distribution  (Ye)  for  five  test 
cases  (i.e.,  damage  scenarios)  for  the  plate  example. 


«  Ci  C2  C  Me(%)  Fe(%)  Ye(%) 

Test  1 

Target 

0.701  0.592  0.511  0.004 

ROM-20 

ROM-40 

ROM-60 

0.544  0.652  0.481  0.003  2.78  3.45  0.37 

0.552  0.610  0.536  0.005  1.53  1.86  0.01 

0.552  0.607  0.521  0.005  1.45  1.57  0.01 

Test  2 

Target 

0.416  0.841  0.832  0.002 

ROM-20 

ROM-40 

ROM-60 

0.310  0.763  0.782  0.003  0.72  3.12  0.04 

0.384  0.803  0.791  0.002  0.52  2.73  0.02 

0.405  0.824  0.801  0.002  0.49  1.67  0.07 

Test  3 

Target 

0.540  0.869  0.264  0.003 

ROM-20 

ROM-40 

ROM-60 

0.407  0.781  0.396  0.005  5.61  4.82  0.32 

0.601  0.855  0.202  0.001  3.52  2.67  0.13 

0.451  0.860  0.236  0.003  1.46  2.22  0.08 

Test  4 

Target 

0.639  0.544  0.647  0.005 

ROM-20 

ROM-40 

ROM-60 

0.558  0.598  0.607  0.003  2.89  5.14  0.48 

0.583  0.576  0.683  0.004  2.4  2.39  0.27 

0.608  0.532  0.651  0.004  2.11  1.5  0.24 

Test  5 

Target 

0.066  0.404  0.448  0.007 

ROM-20 

ROM-40 

ROM-60 

0.048  0.488  0.337  0.005  0.65  4.93  0.07 

0.071  0.381  0.411  0.005  0.48  3.01  0.03 

0.061  0.414  0.491  0.006  0.46  2.47  0.03 

(a)  (b) 

Figure  7:  Spatial  distribution  of  the  Young’s  modulus  from  (a)  the  target  (simulated  experiment) 
and  (b)  the  inverse  characterization  estimate  with  the  ROM  built  from  60  adaptively  generated 
snapshots  for  the  fourth  test  scenario  for  the  plate  example. 


Inverse  Material  Characterization  Using  Gappy  POD  with  Direct  Inversion: 

The  Gappy  POD  process  starts  by  following  the  standard  POD  procedure  to  obtain  a  set 
of  orthogonal  modes  from  a  given  set  of  snapshots.  The  point  in  which  Gappy  POD  diverges 
from  standard  POD  is  how  those  modes  are  utilized.  If  the  full  distribution  (specifically,  full 
spatial  distribution  for  the  cases  herein)  of  a  field  is  available,  the  modal  coefficients  ( a, )  needed  to 
reconstruct  that  field  with  the  POD  modes  can  be  easily  obtained  by  projecting  the  modes  onto 
the  field  as: 


a*  =  /  u  J  •  &  J  d£.  (22) 

j  n  7 

However,  projection  is  no  longer  applicable  to  determine  the  values  of  the  modal  coefficients  to 
reconstruct  the  field  if  the  entire  spatial  distribution  is  not  available.  Thus,  the  objective  of 
Gappy  POD  is  to  provide  a  means  to  reconstruct  the  full  spatial  distribution  of  a  field  using  the 
POD  modes,  but  with  only  a  partial  spatial  distribution  of  the  field  given.  Defining  u  (x)  as  the 
given  partial  distribution  of  the  field  of  interest  such  that  u  (x)  is  (incorrectly)  0  anywhere  data  is 
unavailable,  then  u  (x)  can  be  expressed  in  terms  of  the  corresponding,  but  unknown,  full  spatial 
distribution  as: 


u(x)  —  /3  (x,  u)  u  (x) , 


(23) 


where  j3  (x,  u)  is  a  mask  function  that  is  defined  as  0  where  data  is  unavailable  and  1  where  data  is 
available.  Assuming  that  the  full  spatial  distribution  can  be  approximated  with  the  POD  modes, 
an  approximation  of  u  (x)  can  be  written  in  terms  of  the  POD  modes  as: 


mn 

u*  (x)  =  (3  (x,  u)  (x) .  (24) 

2=1 

Then,  based  upon  a  least-squares  criteria,  the  optimal  set  of  modal  coefficients  to  reconstruct  the 
full  spatial  distribution  of  the  field  can  be  defined  as  that  which  minimizes  an  error  function  of  the 


form: 


2 


e  = 


/ 

Jn 


/3  (x,  u)  u(x)  —  /3  ( x ,  u )  E  ai<j>i(x 

2=1 


dx. 


(25) 


Lastly,  applying  the  necessary  condition  for  extrema  of  a  function  by  setting  the  derivative  of  the 
error  function  with  respect  to  the  modal  coefficients  to  zero,  the  optimal  set  of  modal  coefficients, 
{a},  to  reconstruct  the  full  spatial  distribution  of  the  held  can  be  determined  from  the  solution  of: 


|M]{a}  =  {/}, 


(26) 


where 


and 


Mij 


P  (x,  u)  <f>i(x)  ■  pj(x)  dx. 


fi=/3  (x,  u)  u(x)  •  4>i(x)  dx. 

•/>> 


(27) 

(28) 


Although  potentially  applicable  to  a  variety  of  different  physical  systems,  the  application  of  the 
present  work  is  characterization  of  the  elastic  modulus  distribution  of  a  solid  from  displacement 
measurements  (full-held  displacement  response  once  Gappy  POD  has  been  utilized).  Furthermore, 
the  following  formulation  is  presented  with  respect  to  a  steady-state  dynamic  testing  procedure  (as 
could  be  applicable  to  frequency  response  function-based  evaluation) ,  but  could  easily  be  converted 
to  a  static  problem  by  simply  setting  the  excitation  frequency  to  zero.  Therefore,  the  forward  prob¬ 
lem  and  weak  form  dehned  above  in  Eqns.  (12)  and  (13)  are  still  applicable.  Applying  the  standard 
hnite  element  procedure  and  converting  to  Voigt  notation,  the  domain  can  be  discretized  into  a 
collection  of  elements,  and  the  displacement  and  virtual  displacement  helds  and  their  corresponding 
strain  vectors  can  be  approximated  as: 


u(x,u)  «  [AG(f)]{fT(a;)},  (29) 

Su(x)  ~  [Als(a:)]{5'Ue},  (30) 

(e(f,w)}  «  [£tf(£)]{if  (w)},  (31) 

and 

(fe(f)}  «  [Bz(x)\{6iF},  (32) 

where  [N^x)]  is  the  standard  matrix  of  shape  functions  for  displacement  interpolation  and  [B$  (5)] 
is  the  matrix  of  shape  function  spatial  derivatives  for  strain  interpolation.  Substituting  these 
held  approximations  into  Eqn.  (13),  eliminating  the  arbitrary  virtual  response  held  vector,  and 
assembling  individual  element  contributions,  the  hnal  hnite  element  equations  are  depicted  as: 

[AIM  -  MM  =  {A},  (33) 


where 

[K]=  E  [  [Bu(x)]T[D][Ba(x)}  dx, 

element 


(34) 


(35) 


[M]  =  V  /  p(:r)ta2[A'a(:r)]7'[A'a(:r)]  <ir, 

element  ^ 

{P}=  y;  [  {Nu(x)]tT(x,uj)  dx,  (36) 

Jye 

element  f 

and  [D]  is  elasticity  matrix,  such  that: 

{a(x,ix)}  =  [D]{e(x,u)}  (37) 


The  summation  over  elements  refers  to  the  assembly  process. 

With  the  objective  of  the  inverse  problem  being  characterization  of  the  elastic  modulus  distri¬ 
bution  provided  with  the  entire  displacement  field  everywhere  in  the  domain,  the  first  step  in  the 
inverse  solution  formulation  is  to  separate  the  elastic  modulus  (E(x))  from  the  elasticity  matrix 
as: 

[D]  =  [i Dj\E(x ),  (38) 

where  Di  is  now  only  a  function  of  Poisson’s  ratio  (y).  Applying  the  same  general  weak  form 
procedure  as  was  done  previously  for  displacement,  but  now  for  the  elastic  modulus,  the  inverse 
problem  weak  form  for  the  steady- state  dynamic  boundary  value  problem  can  be  written  as: 


V8E(x)  :  ctj(x,uj)E(x)  dx=  p(x)uj28E(x)  ■  u(x,co)  dx+ 


8E(x)  ■  T(x,  oj)  dx, 


(39) 


where 

{<t/(z»}  =  [D/]{e(x,  w)},  (40) 

and  SE  is  the  virtual  elastic  modulus  vector  (matching  the  dimension  of  the  displacement  field, 
and  therefore,  the  number  of  equilibrium  equations,  even  though  the  modulus  itself  is  a  scalar). 
Now  discretizing  the  domain  into  finite  elements  to  represent  the  elastic  modulus  and  again  using 
Voigt  notation  where  applicable,  the  elastic  modulus  and  virtual  elastic  modulus  vector  and  their 
corresponding  gradients  can  be  approximated  as: 


E(x)  »  [Ayf)])^}. 

(41) 

SE(x)  »  [AT^(f)]{5E'}. 

(42) 

{VE(f)}  «  [B£(f)]{B'>. 

(43) 

{V5E(x)}  «  {Bsg(x)]{SE*}, 

(44) 

where  [Ne(x)]  is  now  the  matrix  of  shape  functions  for  elastic  modulus  interpolation,  [Nsg(x)] 
is  the  expanded  version  (to  match  the  dimensions  of  the  displacement)  of  the  matrix  of  shape 
functions  for  elastic  modulus  interpolation,  and  [Be{x)\  and  [Bsg(x) ]  are  the  respective  matrices 


of  these  shape  function  spatial  derivatives.  Substituting  these  field  approximations  as  well  as  the 
previously-defined  discretization  of  the  given  displacement  field  into  Eqn.  (39),  eliminating  the 
arbitrary  virtual  elastic  modulus  field  vector,  and  assembling  individual  element  contributions,  the 
final  finite  element  equations  for  the  direct  inversion  elastography  problem  are  depicted  as: 

[K,}{E]  =  {/>,}  + (45) 

where 

[*;]=  E  [  iBsdx)}TmBdimwE(t))  dx,  <46) 

element 

[Mi]=  Y,  lp{Z)u2[N&2{x)]T[Nz{x)}dx,  (47) 

element  ^ 

and 

{P,}=  E  /  Wst(Z)]Tf(*,u,)  dx.  (48) 

element  ^7* 

Since  [Kj]  is  typically  non-square  and  Eqn.  (45)  is  typically  an  overdetermined  system  ([Kj]  has 
dimensions  3N  x  N,  where  N  is  the  number  of  nodes  in  the  mesh  if  the  same  mesh  is  used  for  both 
fields),  the  elastic  modulus  cannot  be  estimated  by  simply  inverting  [Kj\.  Thus,  as  is  common, 
a  least-squares  approach  was  used  here  to  solve  Eqn.  (45)  for  {E}.  As  such,  the  nodal  values  of 
elastic  modulus  can  be  determined  as: 

{£}  =  ({K.FIK,})-'  {K,)t  ({P,}  +  \M,\{u}) .  (49) 

The  overall  algorithm  for  direct  inversion  of  a  material  property  distribution  from  partial-field 
response  measurements  with  Gappy  POD  can  be  summarized  as  follows: 

Given:  The  geometry  of  the  structure  of  interest,  the  boundary  conditions  and  partial-field  re¬ 
sponse  measurements  from  a  nondestructive  testing  procedure,  and  any  available  material 
properties. 

Find:  The  unknown  material  property  distribution. 

Step  1  :  Generate  (e.g.,  randomly  or  through  some  other  sampling  procedure)  a  set  of  potential 
distributions  for  the  unknown  material  property,  using  any  information  available  a  priori 
relating  to  the  likely  nature  of  the  unknown  distribution,  and  use  a  forward  analysis  procedure 
to  produce  the  corresponding  full-field  structural  responses  for  each  property  distribution 
from  the  nondestructive  testing  conditions. 

Step  2  :  Calculate  the  POD  modes  from  the  set  of  full-field  structural  responses,  and  select  the 
modes  (based  on  a  user-defined  criteria,  such  as  the  eigenvalue  energy)  to  be  retained  for 
Gappy  POD  field  reconstruction. 

Step  3  :  Reconstruct  the  full-field  structural  response  from  the  given  partial-field  measurements 
with  Gappy  POD. 

Step  4:  Calculate  an  estimate  to  the  unknown  material  property  distribution  using  the  direct 
inversion  procedure  with  the  reconstructed  full-field  structural  response  and  nondestructive 
test  boundary  conditions. 
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Figure  8:  Schematic  for  the  numerically  simulated  examples  representing  characterization  of  an 
elastic  modulus  distribution  with  an  inclusion. 


Examples  were  considered  and  evaluated  to  examine  the  potential  benefits  and  capabilities  of 
using  the  direct  inversion  approach  with  Gappy  POD  to  characterize  the  elastic  modulus  distri¬ 
bution  in  solids  from  partial-field  measurements.  Numerically  simulated  experiments  were  based 
upon  characterization  of  elastic  modulus  distributions  with  circular  inclusions  (hard  or  soft),  as 
shown  schematically  in  Figure  8.  Again,  a  Gaussian  radial  basis  function  representation  was  cho¬ 
sen  to  define  the  localized  elastic  modulus  variations,  as  defined  in  Eqn.  (18).  For  the  inverse 
characterization  process,  it  was  assumed  to  be  known  a  priori  that  the  variation  in  elastic  modulus 
was  similarly  localized  in  nature  (of  course,  with  the  size,  amplitude,  and  location  of  this  variation 
to  be  estimated  by  characterizing  the  entire  spatial  distribution  of  the  elastic  modulus  with  the 
direct  inversion  procedure).  Therefore,  the  process  to  create  the  POD  modes  used  snapshots  gen¬ 
erated  with  this  same  RBF  parameterization  of  the  elastic  modulus  distribution.  For  one  example, 
a  soft  material  block  was  modeled  as  a  50  mm  x  50  mm  square  domain  with  the  bottom  fixed 
to  a  rigid  support  that  was  expected  to  have  a  hard  inclusion.  The  entire  material  (matrix  and 
inclusion)  was  assumed  to  be  known  to  be  nearly  incompressible,  and  a  Poisson’s  ratio  of  0.49999 
was  assigned.  A  simulated  static  test  was  implemented  by  applying  a  0.2  N/mm  (factoring  out  the 
arbitrary  thickness)  excitation  uniformly  to  the  top  surface  of  the  block.  Then,  the  static  vertical 
displacement  response  to  the  loading  was  measured  at  100  uniformly  spaced  discrete  locations,  as 
shown  in  Figure  9.  5%  Gaussian  white  noise  was  added  to  the  measurements  for  this  first  example, 
which  was  deemed  to  be  reasonable  level  of  noise  that  could  be  expected  from  similar  tests  (note 
that  this  level  of  noise  is  commensurate  with  the  highest  levels  of  noise  used  in  prior  referenced 
works  on  direct  inversion  strategies  with  full-field  response  measurements). 

For  the  process  of  generating  the  snapshots  for  POD,  the  elastic  modulus  of  the  background 
material  (i.e.,  matrix  material)  was  assumed  to  be  fixed  at  15  kPa.  Alternatively,  the  parameters 
defining  the  inclusion  based  on  the  RBF  description  were  assumed  to  be  variable.  The  specific 
parameter  values  used  to  create  the  snapshots  were  chosen  arbitrarily  by  uniformly  sampling  the 
space  of  the  four  variable  parameters  (the  two  spatial  coordinates,  amplitude,  and  breadth).  Three 
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Figure  9:  Schematic  of  the  vertical  displacement  sensor  locations  (red  dots)  for  the  example  soft 
matrix  with  a  hard  inclusion. 


Inclusion  Centers 

#  *  * 

#  *  * 

#  *  * 

/ 

\ 

0.3h 

/ 

/ 

\ 

0.2h 

/ 

/ 

\ 

0.2h 

/ 

/ 

\ 

\ 

0.3h 

f 

/  \ 

/  \ 

f'  s 

0.3b 

0.2b 

0.2b 

0.3b 

Figure  10:  Schematic  of  the  nine  inclusion  centers  used  separately  to  generate  the  snapshots  for 
POD  for  the  numerically  simulated  examples. 


values  were  chosen  for  each  spatial  coordinate  of  the  inclusion  center  and  two  values  were  chosen 
each  for  the  amplitude  and  breadth  of  the  inclusion,  and  one  last  scenario  with  no  inclusion  (i.e., 
homogeneous  matrix  material)  was  added,  for  a  total  of  37  parameter  combinations  used  to  create 
snapshots.  Figure  10  shows  the  nine  location  combinations  of  the  inclusion  center  used  to  generate 
the  snapshots.  The  values  of  the  other  two  parameters  used  to  create  the  parameter  combinations 
were  chosen  based  on  an  expectation  of  what  the  lower  and  upper-end  would  be  for  the  application, 
using  1  and  3  for  the  amplitude  parameter  (i.e.,  modulus  at  inclusion  center  of  30  kPa  and  60  kPa) 
and  5  mm  and  15  mm  for  the  breadth  parameter. 

Figure  11  shows  a  representative  example  trial  of  a  simulated  displacement  response  field  from  a 
randomly  generated  parameter  set  and  including  the  5%  Gaussian  white  noise  in  comparison  to  the 
displacement  response  field  reconstructed  from  the  100  response  measurements  of  that  response 
field  with  the  Gappy  POD  procedure.  In  general,  the  Gappy  POD  reconstruction  of  the  dis¬ 
placement  response  from  partial-field  measurements  was  found  to  be  accurate,  producing  response 
distributions  that  were  nearly  identical  to  the  full  simulated  responses,  with  errors  consistent  with 
the  example  shown,  which  had  relative  L2  and  errors  in  the  displacement  reconstruction  of  7.4% 


,-4 


(a) 


,-4 


3.00x10 

1.50x10 

0.00 

-1.50x10" 

-3.00x10" 


3.00x10" 


(c) 


8.00xl0-4 
6.00x1  O'4 
4.00x1  O'4 
2.00x1  O'4 
0.00 


8.00xl0-4 
6.00x1  O'4 
4.00x1  O'4 
2.00x1  O'4 
0.00 


(d) 


Figure  11:  Representative  example  of  the  (a)  horizontal  and  (b)  vertical  components  of  a  simulated 
experimental  displacement  field  with  5%  Gaussian  white  noise  and  the  (c)  horizontal  and  (d) 
vertical  components  of  the  corresponding  reconstructed  displacement  field  from  Gappy  POD  with 
only  the  discrete  measurement  data  (color  contours  in  units  of  m)  for  the  example  soft  matrix  with 
a  hard  inclusion. 
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Figure  12:  Representative  example  of  (a)  the  elastic  modulus  distribution  used  to  simulate  ex¬ 
perimental  measurements  (i.e.,  the  target  modulus  distribution)  and  (b)  the  corresponding  elastic 
modulus  distribution  estimated  with  the  direct  inversion  approach  with  Gappy  POD  full-field  dis¬ 
placement  reconstruction  from  the  partial-field  measurements  (color  contours  in  units  of  Pa)  for 
the  example  soft  matrix  with  a  hard  inclusion. 


and  18.5%,  respectively.  More  importantly,  Figure  12  shows  a  representative  example  of  an  elastic 
modulus  distribution  used  to  generate  simulated  experimental  measurements  (i.e.,  target  modulus 
distribution)  and  the  corresponding  elastic  modulus  distribution  estimated  by  the  direct  inversion 
procedure  with  Gappy  POD  displacement  reconstruction.  There  was  a  noticeable  amount  of  er¬ 
ror  in  the  modulus  reconstruction,  and  the  relative  L2  and  L~^  errors  in  the  modulus  estimation 
compared  to  the  target  distribution  were  21%  and  43%,  respectively.  However,  the  localization 
of  the  modulus  distribution  was  accurate,  indicating  a  single  harder  region  within  the  solid,  and 
the  matrix  modulus  value  and  the  maximum  inclusion  modulus  value  were  nearly  identical  to  the 
target  distribution.  Again,  it  should  be  noted  that  the  direct  inversion  process  does  not  restrict  the 
distribution  of  the  modulus  to  be  localized,  which  emphasizes  the  significance  of  having  recovered 
the  correct  localization  with  the  inversion  procedure.  Moreover,  the  accuracy  of  the  modulus  recon¬ 
struction  was  commensurate,  if  not  qualitatively  better,  than  in  alternate  works  in  the  literature  on 
direct  inversion  techniques  [10,  11],  with  those  approaches  using  full-field  response  measurements, 
rather  than  the  partial-field  measurements  used  here.  To  further  expand  on  the  benefits  of  the 
direct  inversion  with  Gappy  POD,  Figure  13  shows  the  elastic  modulus  distribution  estimated  by 
the  direct  inversion  procedure  with  the  original  simulated  full-field  displacement  with  the  added 
5%  Gaussian  white  noise.  Clearly,  the  direct  inversion  procedure  applied  to  the  noisy  data  was 
unable  to  remotely  come  close  to  estimating  the  target  elastic  modulus  distribution,  indicating  the 
significance  of  the  noise  filtering  capability  of  the  Gappy  POD  field  reconstruction  prior  to  direct 
inversion.  For  all  test  cases  investigated  the  Gappy  POD  with  direct  inversion  procedure  was  able 
to  reconstruct  the  full-field  displacement  and  then  estimate  the  elastic  modulus  distribution  with 
an  accuracy  consistent  with  the  representative  example  shown.  In  particular,  the  capability  of 
the  Gappy  POD  with  direct  inversion  procedure  to  correctly  localize  and  accurately  estimate  the 
magnitude  of  the  elastic  modulus  distribution  remained  consistent. 
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Figure  13:  Example  of  an  elastic  modulus  distribution  estimated  with  the  direct  inversion  approach 
applied  directly  to  the  full-field  simulated  experimental  data  with  noise  (i.e. ,  without  using  Gappy 
POD)  (color  contours  in  units  of  Pa)  for  the  example  soft  matrix  with  a  hard  inclusion. 


Multi- Objective  Optimization  for  Computationally  Efficient  Nondestructive  Evalu¬ 
ation: 

This  portion  of  the  work  also  focused  on  NDE  of  solid  continua  and  utilized  the  common  ap¬ 
proach  of  casting  the  inverse  characterization  problem  as  an  optimization  problem  to  determine 
the  unknown  parameters  of  the  structure  to  minimize  the  difference  between  the  measured  re¬ 
sponse  and  those  predicted  by  a  numerical  representation  of  the  structure  (e.g.,  finite  element  or 
boundary  element  analysis)  subject  to  the  nondestructive  testing  conditions.  This  type  of  opti¬ 
mization  strategy  to  approximate  the  solution  of  inverse  characterization  problems  provides  several 
benefits,  including  the  capability  to  address  a  variety  of  physical  processes/measurements,  even  si¬ 
multaneously,  and  providing  quantitative  results.  However,  there  are  also  inherent  challenges  in 
addition  to  the  ubiquitous  ill-posedness  associated  with  inverse  problems,  most  prevalent  being 
the  large  computational  expense  of  this  optimization  approach.  The  key  to  using  multi-objective 
optimization  for  these  inverse  problems  is  as  simple  as  dividing  up  the  measurement  components 
or  spatial  distribution  into  several  separate  objective  functionals  to  be  minimized  simultaneously, 
but  separately,  which  could  be  viewed  as  (e.g.,  for  the  case  of  dividing  the  spatial  distribution  of 
the  measurements): 


Minimize  < 


\\Rsim(a,x) 
||  Rsim(a,x) 


a 


Rexp{x)  ||ri 
Rexp{x)  ||r2 


(50) 


J\Rsim(a,x)  -  Rexp(x)\\Tn 


where  a  is  the  vector  of  parameters  to  be  determined  to  characterize  the  desired  structural  proper¬ 
ties,  Rsim  is  the  simulated  response  field  to  estimate  the  NDT  response  for  a  given  set  of  parameters, 
Rexp  is  the  experimentally  measured  response  field  (i.e.,  optimization  target),  ||  •  ||r  is  some  suitable 
metric  norm  with  respect  to  the  domain  of  the  nondestructive  measurements,  r.t  is  the  ith  subdi¬ 
vision  of  the  domain  of  the  response  field  measurements  obtained  from  nondestructive  testing  and 
n  is  the  total  number  of  subdivisions.  Then,  any  preferred  multi-objective  optimization  algorithm 
can  be  employed  to  determine  the  Pareto  front  for  Eqn.  (50),  which  can  be  thought  of  as  the  set  of 


all  possible  solutions  to  the  inverse  problem  that  have  a  lower  value  for  at  least  one  of  the  separate 
objective  functionals  in  comparison  to  any  other  solution  estimate  seen  throughout  the  optimiza¬ 
tion  process.  This  work  employed  a  controlled  elitist  multi-objective  genetic  algorithm  (CEMGA) 
[12,  13]  to  determine  the  Pareto  front  in  the  example  cases  considered  therein.  When  the  multi¬ 
objective  optimization  process  is  complete,  a  single  solution  estimate  for  the  inverse  problem  can 
be  attained  through  some  chosen  final  decision  criteria,  such  as  the  minimum  sum  of  all  objective 
functionals,  such  as: 

n 

Minimize  II Rsim(a,  x)  -  Rexp(x )  ||r, ,  (51) 

<3e{<5i}i=i  .=1 

where  {cq}f=1  is  the  set  of  p  potential  solutions  identified  as  part  of  the  Pareto  front.  Apply¬ 
ing  multi-objective  optimization  rather  than  the  standard  single  objective  optimization  to  inverse 
characterization  problems  was  found  to  be  a  simple  yet  effective  approach  to  substantially  re¬ 
duce  the  computational  expense  and  improve  the  consistency  of  the  solution  accuracy  for  inverse 
characterization. 

The  primary  feature  of  the  multi-objective  optimization  that  leads  to  improved  inverse  solution 
capabilities  is  that  substantial  diversity  of  the  solution  estimates  is  maintained  throughout  the 
search  process  by  evolving  a  set  of  optima  (i.e.,  the  Pareto  front)  rather  than  a  single  optimum 
throughout  an  iterative  optimization  process.  By  maintaining  diversity,  the  multi-objective  opti¬ 
mization  process  is  uniquely  able  to  traverse  the  large  parameter  search  spaces  that  are  typical 
of  inverse  characterization  problems  efficiently  and  consistently,  avoiding  stalling  and  convergence 
to  local  minima.  An  additional  benefit  of  the  diversity  in  the  solution  estimates  provided  by 
multi-objective  optimization  is  the  resulting  improvement  in  the  ability  to  reveal  the  variety  of 
solutions  that  may  exist  for  ill-posed  (particularly  non- unique)  problems.  As  a  direct  reason  for 
non-uniqueness  can  be  insufficiency  of  the  parameterization  of  the  properties  to  be  determined, 
the  solution  diversity  provided  by  multi-objective  optimization  can  thus  be  assumed  to  be  able  to 
provide  insight  into  the  changes  to  the  parameterization  necessary  to  subsequently  produce  more 
unique  and  accurate  inverse  solutions.  To  provide  additional  context,  the  following  discussion 
will  provide  example  scenarios  based  upon  the  class  of  NDE  problems  relating  to  characterizing 
an  unknown  quantity  of  localized  changes  in  properties  (e.g.,  damage  or  defect  characterization). 
However,  this  approach  should  be  able  to  be  similarly  implemented  for  a  variety  of  inverse  char¬ 
acterization  problems,  particularly  those  for  which  some  property  of  the  unknown  field  is  known  a 
priori  to  be  (semi-)  localized  in  the  parameter  space. 

The  overall  structure  of  the  optimization-type  NDE  algorithm  incorporating  an  adaptive  self- 
evolving  parameterization  approach  with  multi-objective  optimization  is  shown  in  Figure  14.  The 
core  hypothesis  of  the  self-evolving  parameterization  component  of  the  algorithm  developed  is  that 
the  distribution  of  solutions  in  the  parameter  space  produced  through  multi-objective  optimization 
provides  guidance  as  to  whether  a  specific  physical  parameter  should  be  expanded.  In  the  localized 
property  characterization  context,  in  which  the  primary  method  to  expand  the  parameterization 
could  be  to  increase  the  number  (n)  of  basis  functions  with  compact  (or  semi-compact)  support 
used  to  define  the  property  distribution  along  with  their  associated  unknown  parameters  to  be 
determined  by  the  inversion,  the  self-evolving  parameterization  component  of  the  characterization 
algorithm  could  be  implemented  as  follows: 


Figure  14:  Flowchart  of  the  multi-objective  optimization-type  NDE  algorithm  with  adaptive  self- 
evolving  parameterization. 


Given  -  The  Pareto  front  of  potential  solutions  to  the  inverse  characterization  problem  (e.g., 
spatial  coordinates  of  the  centroid  for  each  localized  property  change  and  any  associated  pa¬ 
rameters)  subject  to  the  current  value  of  n  (i.e. ,  the  number  of  compactly  or  quasi-compactly 
supported  basis  functions  used  to  define  the  localized  change). 

Step  1  -  Identify  the  n  +  1  regions  of  localized  property  change  whose  centroids  are  separated 
by  the  largest  Euclidean  distance  from  the  entire  Pareto  set  of  solutions  (noting  that  each 
solution  set  could  contain  multiple  localized  property  changes  depending  on  the  value  of  n), 
referred  to  as  the  n  +  1  “Parameterization  Poles” . 

Step  2  -  Identify  and  average  all  regions  of  localized  property  changes  that  overlap  with  each 
Parameterization  Pole  to  produce  the  n+1  “Cluster  Means”  of  the  localized  property  changes. 
Step  3  -  Do  any  of  the  Cluster  Means  overlap? 

Yes  — >■  STOP  (do  not  update  the  parameterization  further). 

No  — >  SET  n  =  n  +  1  and  GO  TO  Step  1. 

To  examine  the  capability  of  the  self-evolving  parameterization  approach  utilizing  multi-objective 
optimization  for  NDE  to  efficiently  and  accurately  characterize  localized  property  changes  in  solid 
continua  several  simulated  examples  of  damage  characterization  within  structural  steel  plates  (as 
could  potentially  be  affected  by  erosion)  were  considered.  More  specifically,  the  example  NDE 
cases  sought  to  characterize  the  size  and  location  of  circular  regions  of  material  loss  within  the 
steel  plates  considered.  Figure  15  shows  a  schematic  for  one  set  of  test  cases,  which  consisted  of  an 
arbitrarily  thin  lmxlm  square  section  with  the  bottom  fixed  to  a  rigid  support.  The  simulated 
NDT  consisted  of  applying  a  1  kN/m  (factoring  out  the  arbitrary  thickness)  harmonic  excitation 


Figure  15:  Schematic  of  the  damaged  square  plate  example. 


to  the  top  surface  of  the  plate  at  an  excitation  frequency  of  20  Hz,  and  then  measuring  the  vertical 
and  horizontal  displacements  at  99  equally-spaced  increments  along  the  left  and  right  surfaces. 
For  both  generating  the  experimental  data  and  simulating  the  forward  problem  during  the  inverse 
solution  process  the  structures  were  assumed  to  behave  linearly  and  be  defined  by  steady-state 
dynamic  plane  stress  solid  mechanics,  and  all  analyses  were  performed  using  the  finite  element 
method.  The  inverse  problems  to  determine  the  parameters  defining  the  damage  in  the  example 
structures  were  cast  in  the  form  of  the  following  multi-objective  optimization  problem,  arbitrarily 
having  selected  to  divide  the  displacement  measurements  into  four  groupings: 
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where  a  is  the  vector  containing  the  parameters  of  the  unknown  damage  to  be  determined  in  the 
inverse  problem,  U'{'-p .  U^p ,  U^p  and  U'4JP  are  the  experimentally  measured  displacement  at  the 
jth  measurement  location  in  the  four  groupings,  and  t/j’j'" ,  ,  U4"n  and  U4"n  are  the  numerically 

simulated  displacement  at  the  jth  measurement  location  in  the  four  groupings.  The  four  objective 
functions  for  this  example  were  simply  defined  by  dividing  the  measurements  with  respect  to  the 
two  sides  and  the  two  directional  components. 

One  of  the  simulated  scenarios  to  test  the  capabilities  of  the  NDE  algorithm  with  self-evolving 
parameterization  was  a  case  with  two  actual  damage  regions  in  the  simulated  experiment.  Figure 
16  shows  a  representative  example  from  five  trials  of  the  inverse  solution  process  of  the  Pareto 
front  solution  estimates  and  Figure  17  shows  the  corresponding  measurement  error  for  the  four 
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Figure  16:  Representative  results  of  the  Pareto  front  damage  region  solution  estimates  obtained 
from  the  NDE  algorithm  with  self-evolving  parameterization  after  (a)  the  first  iteration  with  a 
single  damage  region  parameterization  and  (b)  the  second  iteration  for  which  the  parameterization 
had  been  automatically  updated  to  two  damage  regions,  along  with  the  best  individual  from  each 
Pareto  front,  and  compared  to  the  actual  (experimental)  two  damage  regions  for  the  example 
square  plate. 


objectives  for  the  “best”  individual  from  each  of  the  Pareto  fronts  at  each  solution  iteration  (as 
the  parameterization  evolved),  and  Figure  18  shows  the  final  solution  estimate  from  the  converged 
algorithm  for  this  example  with  two  actual  damage  regions.  The  final  solution  estimates  for  all 
trials  of  the  test  cases  were  similarly  accurate  as  this  example  shown,  with  both  the  sizes  and  loca¬ 
tions  of  the  damage  regions  being  relatively  accurate.  The  ability  of  the  adaptive  inverse  solution 
process  to  consistently  and  efficiently  determine  the  exact  number  of  damage  regions  and  provide 
relatively  accurate  estimations  of  the  location  and  size  of  these  damages  based  solely  on  surface 
measurements  is  significant.  Moreover,  the  NDE  algorithm  with  self-evolving  parameterization  was 
found  to  be  robust  to  potential  noise  in  the  measurement  data. 
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